#!/usr/bin/perl -w
use Getopt::Long;
use FindBin;
use strict;
use DBI;
use DBD::Pg;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


my $usage = <<'USAGE';

USAGE:

    panda_omim_codes.pl --morbid_omim FILE --genemap_omim FILE
                        --db_omim FILE [OPTIONS]...  DB_NAME

    MANDATORY arguments:

        --morbid_omim FILE   the OMIM morbid map available from
                                   ftp://ftp.ncbi.nih.gov/repository/OMIM/morbidmap
        --genemap_omim FILE  the OMIM gene map available from
                                   ftp://ftp.ncbi.nih.gov/repository/OMIM/genemap
        --omim_db FILE       the full OMIM database available as a compressed archive from
                                   ftp://ftp.ncbi.nih.gov/repository/OMIM/omim.txt.Z

    Optional arguments:

        [--user USER_NAME]   the postgreSQL user name of the owner/panda administrator
        [--host STRING]      the postgreSQL database server
        [--password STRING]  the postgreSQL database password
        [--verbose]          print more messages showing progress


        [--help]


RESULTS:


    This script creates the panda system with DB_NAME as the PostgreSQL database
                behind panda.

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


my $user			= 'postgres';
my $password		= undef;
my $host			= '';
my $verbose;
my $genemap_omim	= '';
my $morbid_omim 	= '';
my $omim_db     	= '';
my $help			= '';
my $dbname			= '';

GetOptions (
			
			'genemap_omim=s'	=> \$genemap_omim,
			'morbid_omim=s'		=> \$morbid_omim,
			'omim_db=s'			=> \$omim_db,
			'help'				=> \$help,
			'user=s'			=> \$user,
			'password'			=> \$password,
			'host=s'			=> \$host,
			'verbose'	 		=> \$verbose,
			);


die $usage if ($help);
$dbname			= $ARGV[0] if ($ARGV[0]);
die "\nError:\n\tMissing database name.\n$usage" unless ($dbname);
die "\nError:\n\tMissing genemap_omim name.\n$usage" unless ($genemap_omim);
die "\nError:\n\tMissing morbid_omim  name.\n$usage" unless ($morbid_omim );
die "\nError:\n\tMissing omim_db      name.\n$usage" unless ($omim_db     );
die "\nError:\n\tMissing user name.\n$usage" unless ($user);
$host		= ';host=' 	. $host if ($host);



#_________________________________________________________________________________________

#	Parse OMIM codes from the gene map

#_________________________________________________________________________________________
	
my %MIM_codes_gm;
print STDERR "    Parsing OMIM codes from the OMIM gene map...\n"
	if ($verbose);

{
	open GENEMAP_OMIM, $genemap_omim or die "Error\n\tCould not open the [$genemap_omim] file\n$@\n";
	my $cnt_MIMs;
	while (<GENEMAP_OMIM>)
	{
		if(!
			/	^							# start at beginning of line
			(?:[^|]*\|){9}				# ignore first 9 fields
			([^|]*)\|					# save MIM code field (not including | terminator)
			(?:[^|]*\|){3}				# ignore next 3 fields
			([^|]*)						# save disease field
							/x)
		{
			next;
		}

		if ($1 ne '' && $2 ne '')
		{
			my ($MIM_code, $dis) = ($1,$2);
			
			if ($dis !~ /^(		[?\[]		|		# not begin with ? or [
								(\{\?)		)	/x)	# or {?
			{
				++$MIM_codes_gm{$MIM_code};
				++$cnt_MIMs;
			}
		}
	}
	if ($verbose)
	{
		print STDERR "\t", scalar keys %MIM_codes_gm,
						"\tMIM codes were retrieved from the gene map.\n";
		print STDERR "\t\t($cnt_MIMs gene map entries were associated with MIM codes.)\n";
	}
}	





#_________________________________________________________________________________________

#	Parse OMIM codes from the morbid map

#_________________________________________________________________________________________
my %MIM_codes_mm;
print STDERR "    Parsing OMIM codes from the morbid map...\n"
	if ($verbose);

{
	open MORBID_OMIM , $morbid_omim  or
				die "Error\n\tCould not open the [$morbid_omim] file\n$@\n";
	my $cnt_MIMs;
	while (<MORBID_OMIM>)
	{
		if(!
		/	^							# start at beginning of line
			(?:[^|]*\|){2}				# ignore first 2 fields
			(\d+)\|					    # save MIM code field (not including | terminator)
							/x)
		{
			next;
		}
		++$MIM_codes_mm{$1};
		++$cnt_MIMs;
	}
	if ($verbose)
	{
		print STDERR "\t", scalar keys %MIM_codes_mm, "\tMIM codes were retrieved from the morbid map.\n";
		print STDERR "\t\t($cnt_MIMs morbid map entries were associated with MIM codes.)\n";
	}
}	



#_________________________________________________________________________________________

#	Parse OMIM codes associated with variants from the OMIM gene map

#_________________________________________________________________________________________
my %MIM_codes_db;
print STDERR "    Parsing OMIM codes associated with variants from the OMIM database...\n"
	if ($verbose);
{
	open DB_OMIM, "lookup_variants_in_omim " . ($verbose ? "--v " : "") . $omim_db . "|"
			or die "Error\n\t".
					"Could not run [lookup_variants_in_omim] to parse OMIM codes ".
					"associated with variants in the OMIM db [$omim_db]\n".
					"lookup_variants_in_imim can be built by running (in this directory):\n".
					"gmake -f /bin/lookup_variants_in_omim/lookup_variants_in_omim.mak CFG=Release\n".
					"[$@]\n";
	my $cnt_MIMs;
	while (<DB_OMIM>)
	{
		chomp;
		++$MIM_codes_db{$_};
		++$cnt_MIMs;
	}
}	







sub uniquify(\@)
{
	my %hash;
	@hash{@{$_[0]}} = ();
	@{$_[0]} = keys %hash;	
	return @{$_[0]};
}	

my %MIM_codes;
$MIM_codes{$_}++ for (keys %MIM_codes_mm);
$MIM_codes{$_}++ for (keys %MIM_codes_db);
$MIM_codes{$_}++ for (keys %MIM_codes_gm);
if ($verbose)
{
	print STDERR "    Unique OMIM codes associated with disease:\n\t", ('_' x 5), "\n\t",
					scalar keys %MIM_codes, "\tin total.\n";
}


# connect to database
my $dbh = DBI->connect( "dbi:Pg:dbname=$dbname $host",
                        $user,
                        $password,
						{
                          RaiseError => 0,
                          AutoCommit => 1,
						  PrintError => 0
                        }
                      )
		or die "Database connection not made: $DBI::errstr";

# prepare SQL
$dbh->do("Create table omim_morbid_codes( code INTEGER);");
my $sth = $dbh->prepare("COPY newtaxid_from_taxname ( tax_name_hash, tax_id) FROM STDIN");
$sth->execute();

print STDERR "\t[";
my $cnt_names = 0;
while (<NAMES>)
{
	# put name into panda
	next unless /^(\d+)\t\|\t([^\t]+)\t\|/;
	if (++$cnt_names % 4000 == 0)
	{
		print STDERR "\n\t" if (($cnt_names / 4000) % 50 == 0);
		print STDERR '.';
	}
	my @digest = unpack("L*", $md4->hash($2));
	my $line = '{'.join(',', @digest)."}\t$1\n";
	$dbh->func($line, 'putline');
	
	# save name for lookup tree later
	#1	|	root	|		|	scientific name	|
	next unless /^(\d+)\t\|\t([^\t]+)\t\|[^\|]+\|\tscientific name/;
	$names{$1} = $2;
}
$dbh->func("\\.\n", 'putline');
$dbh->func('endcopy');
print STDERR "]\n";

$dbh->disconnect();

